home
***
CD-ROM
|
disk
|
FTP
|
other
***
search
/
Info-Mac 4
/
Info_Mac IV CD-ROM (Pacific HiTech Inc.)(August 1994).iso
/
Science
/
RLaB
/
help
/
filter
< prev
next >
Wrap
Text File
|
1994-04-25
|
3KB
|
115 lines
Synopsis: Digital filter implementation
Syntax: filter ( B , A , X )
filter ( B , A , X , Zi )
Description:
Filter is an implementation of the standard difference
equation:
y[n] = b(1)*x[n] + b(2)*x[n-1] + ... b(nb+1)*x[n-nb]
- a(2)*y[n-1] - ... a(na+1)*y[n-na]
The filter is implemented using a method described as a
"Direct Form II Transposed" filter. More for information see
Chapter 6 of "Discrete-Time Signal Processing" by Oppenheim
and Schafer.
The inputs to filter are:
B The numerator coefficients, or zeros of the
system transfer function. The coefficients are
specified in a vector like:
[ b(1) , b(2) , ... b(nb) ]
A The denominator coefficients, or the poles of
the system transfer function. the coefficients
are specified in a vector like:
[ a(1) , a(2) , ... a(na) ]
X A vector of the filter inputs.
Zi [Optional] The initial delays of the filter.
The filter outputs are in a list with element names:
y The filter output. y is a vector of the same
dimension as X.
zf A vector of the final values of the filter
delays.
The A(1) coefficient must be non-zero, as the other
coefficients are divided by A(1).
Below is an implementation of filter() in a r-file - it is
provided for informational usage only.
#--------------------------------------------------------------
# Simplistic version of RLaB's builtin function filter()
# Y = filter ( b, a, x )
# Y = filter ( b, a, x, zi )
#
rfilter = function ( _b , _a , x , zi )
{
local (M, N, NN, a, b, k, n, ntotal, v, vi, y)
# Copy _a and _b cause we are going to modify them
a = _a;
b = _b;
ntotal = x.nr * x.nc;
M = b.nr * b.nc;
N = a.nr * a.nc;
NN = max ([ M, N ]);
y = zeros (x.nr, x.nc);
# Fix up pole and zero vectors.
# Make them the same length, this makes
# filter's job much easier.
if (N < NN) { a[NN] = 0; }
if (M < NN) { b[NN] = 0; }
# Adjust filter coefficients
if (a[1] == 0) { error ("rfilter: 1st A term must be non-zero"); }
a[2:NN] = a[2:NN] ./ a[1];
b = b ./ a[1];
# Create delay vectors and load inital delays.
# Add an extra term to vi[] to make filter's
# job a little easier. This extra term will
# always be zero.
v = zeros (NN, 1);
vi = zeros (NN+1, 1);
if (exist (zi))
{
vi[1:NN] = zi;
}
#
# Do the work...
#
for (n in 1:ntotal)
{
v[1] = b[1]*x[n] + vi[2];
y[n] = v[1];
for (k in 2:NN)
{
v[k] = b[k]*x[n] - a[k]*v[1] + vi[k+1];
vi[k] = v[k];
}
}
return << y = y; zf = v >>;
};
#--------------------------------------------------------------